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LATTICE  BOLTZMANN  METHOD  FOR  3-D  FLOWS  WITH  CURVED  BOUNDARY 


RENWEI  MEI*,  WEI  SHYY+,  DAZIII  YU*,  AND  LI-SHI  LUO§ 

Abstract.  In  this  work,  we  investigate  two  issues  that  are  important  to  computational  efficiency  and 
reliability  in  fluid  dynamics  applications  of  the  lattice  Boltzmann  equation  (LBE):  (1)  Computational  stabil¬ 
ity  and  accuracy  of  different  lattice  Boltzmann  models  and  (2)  the  treatment  of  the  boundary  conditions  on 
curved  solid  boundaries  and  their  3-D  implementations.  Three  athermal  3-D  LBE  models  (D3Q15,  D3Q19, 
and  D3Q27)  are  studied  and  compared  in  terms  of  efficiency,  accuracy,  and  robustness.  The  boundary  treat¬ 
ment  recently  developed  by  Filippova  and  Hanel  and  Mei  et  al.  in  2-D  is  extended  to  and  implemented  for 
3-D.  The  convergence,  stability,  and  computational  efficiency  of  the  3-D  LBE  models  with  the  boundary 
treatment  for  curved  boundaries  were  tested  in  simulations  of  four  3-D  flows:  (1)  Fully  developed  flows  in 
a  square  duct,  (2)  flow  in  a  3-D  lid-driven  cavity,  (3)  fully  developed  flows  in  a  circular  pipe,  and  (4)  a 
uniform  flow  over  a  sphere.  We  found  that  while  the  fifteen-velocity  3-D  (D3Q15)  model  is  more  prone  to 
numerical  instability  and  the  D3Q27  is  more  computationally  intensive,  the  D3Q19  model  provides  a  balance 
between  computational  reliability  and  efficiency.  Through  numerical  simulations,  we  demonstrated  that  the 
boundary  treatment  for  3-D  arbitrary  curved  geometry  has  second-order  accuracy  and  possesses  satisfactory 
stability  characteristics. 

Key  words,  lattice  Boltzmann  method,  boundary  condition  for  curved  geometries,  accuracy,  3-D  flows, 
Navier-Stokes  equations 

Subject  classification.  Fluid  Dynamics 

1.  Introduction. 

1.1.  Basic  Notion  of  the  Lattice  Boltzmann  Equation.  In  one  fashion  or  another,  conventional 
methods  of  computational  fluid  dynamics  (CFD)  compute  pertinent  flow  fields,  such  as  velocity  u  and 
pressure  p,  by  numerically  solving  the  Navier-Stokes  equations  in  space  x  and  time  t  [23,  9,  27].  In  contrast, 
various  kinetic  methods  use  the  transport  equation,  or  the  Boltzmann  equation  in  particular,  for  various 
problems  in  fluid  dynamics.  The  Boltzmann  equation  deals  with  the  single  particle  distribution  function 
f(x,  £,  t),  where  £  is  the  particle  velocity,  in  phase  space  {x,  £)  and  time  t.  Recently,  the  method  of  the 
lattice  Boltzmann  equation  (LBE)  [5,  24,  3,  6]  has  become  an  alternative  to  the  conventional  CFD  methods 
employing  Navier-Stokes  equations.  The  theoretical  premises  of  the  LBE  method  are  that  (1)  hydrodynamics 
is  insensitive  to  the  details  of  microscopic  physics,  and  (2)  hydrodynamics  can  be  preserved  so  long  as  the 
conservation  laws  and  associated  symmetries  are  respected  in  the  microscopic  or  mesoscopic  level.  Therefore, 
the  computational  advantages  of  the  LBE  method  are  attained  by  drastically  reducing  the  particle  velocity 
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space  £  to  only  a  very  few  discrete  points  without  seriously  degrading  hydrodynamics.  This  is  possible 
because  the  LBE  method  rigorously  preserves  the  hydrodynamic  moments  of  the  distribution  function  /, 
such  as  mass  density  and  momentum  fluxes,  and  the  necessary  symmetries  [13,  14,  1], 

One  popular  kinetic  model  is  the  Boltzmann  equation  with  the  single  relaxation  time  approximation  [4] : 

(i-i)  at/  +  £-v/  =  -l[/-/(°)], 

where  £  is  the  particle  velocity,  /(0)  is  the  equilibrium  distribution  function  (the  Maxwell-Boltzmann  distri¬ 
bution  function),  and  A  is  the  relaxation  time.  The  mass  density  p  and  momentum  density  pu  are  the  first 
( D  + 1)  hydrodynamic  moments  of  the  distribution  function  /  and  /(0) ,  where  D  is  the  dimension  of  velocity 
space. 

To  solve  for  /  numerically,  Eq.  (1.1)  is  first  discretized  in  the  velocity  space  £  using  a  finite  set  of 
velocities  {£Q}  without  affecting  the  conserved  hydrodynamic  moments  [13,  14,  1], 

(1-2)  dtfa  +  *a-Vf«  =  -j[fa-f!?)]- 

In  the  above  equation,  fa(x,  t)  =  /( x,  £a,t)  and  fa\x,  t )  =  /(0)  (*,£„,  t)  are  the  distribution  function 
and  the  equilibrium  distribution  function  of  the  a-th  discrete  velocity  £a ,  respectively.  The  nine- velocity  (or 
9-bit)  LBE  model  on  the  2-D  square  lattice,  denoted  as  D2Q9  model,  has  been  widely  used  for  simulating  2-D 
flows.  For  3-D  flows,  there  are  several  cubic  lattice  models,  such  as  the  fifteen- velocity  (D3Q15),  nineteen- 
velocity  (D3Q19),  and  twenty-seven-velocity  (D3Q27)  models  [12],  which  have  been  used  in  the  literature. 
All  three  models  have  a  rest  particle  (with  zero  velocity)  in  the  discretized  velocity  set  {£Q}.  A  minor 
variation  of  those  models  is  to  remove  the  rest  particles  from  the  discrete  velocity  set;  the  resulting  models 
are  known  as  the  D3Q14,  D3Q18,  and  D3Q26  models,  respectively.  The  LBE  models  with  a  rest  particle 
generally  have  better  computational  stability.  For  athermal  fluids,  the  equilibrium  distributions  for  D2Q9, 
D3Q15,  D3Q19,  and  D3Q27  models  are  all  of  the  form  [14] 

l  +  ^ea-u  +  ^(ea-u)  -—(u-u)  , 

where  wa  is  a  weighting  factor  and  ea  is  a  discrete  velocity,  c  =  Sx/St  is  the  lattice  speed,  and  Sx  and  St 
are  the  lattice  constant  and  the  time  step,  respectively.  (The  values  of  the  weighting  factor  wa  for  D3Q15, 
D3Q19,  and  D3Q27  models  and  the  diagrams  illustrating  the  lattice  structure  for  D3Q15  and  D3Q19  models 
are  given  in  the  Appendix.)  It  can  be  shown  that  fa9 '  is  in  fact  a  Taylor  series  expansion  of  the  Maxwellian 
/a°'  [13,  14].  This  approximation  of  /q0)  by  the  above  fa  makes  the  method  valid  only  in  the  incompressible 
limit  u/c  — >  0. 

With  the  velocity  space  discretized,  the  hydrodynamic  moments  of  /  and  /(0)  are  evaluated  by  the 
following  quadrature  formulas: 

(!-4a)  P  =  ^fa  =  ^2f{a\ 

a  a 

(1  -4b)  P«  =  £ea/a  =  ]>>/r 

a  a 

The  speed  of  sound  of  the  above  3-D  LBE  models  is  cs  =  c/\/3  and  the  equation  of  state  is  that  of  an  ideal 
gas  p  =  pc2s .  The  viscosity  of  the  fluid  is  v  =  Ac^ . 

Equation  (1.2)  is  often  discretized  in  space  x  and  time  t  into 

(1.5)  fa(Xi  +  eaSt,  t  +  St )  +  fa(Xi,  t )  =  \ja(Xi,  t)  -  f(°\xt,  t)]  , 


(1.3) 


&q)  =  wn 
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where  r  =  A /St.  This  is  the  lattice  Boltzmann  equation  with  the  Bhatnagar-Gross-Krook  (BGK)  approx¬ 
imation  [4]  and  is  often  referred  to  as  the  LBGK  model  [5,  24].  The  viscosity  in  the  NS  equation  derived 
from  Eq.  (1.5)  is 

(1.6)  v  =  (r  —  1/2  )c2s5t. 


This  choice  for  the  viscosity  makes  the  LBGK  scheme  formally  a  second  order  method  for  solving  incom¬ 
pressible  flows  [14].  The  positivity  of  the  viscosity  requires  that  r  >  1/2.  Equation  (1.5)  can  be  solved  in 
the  following  two  steps: 

(1.7a)  collision  step:  fa&u  t)  =  fa{xu  t)  -  ^  [/„(*,;,  t)  -  fi°](Xi,  £)]  , 

(1.7b)  streaming  step:  fa{Xi  +  eaSt,  t  +  St)  =  fa(xi,  t), 


where  fa  denotes  the  post-collision  state  of  the  distribution  function.  It  is  noted  that  the  collision  step  is 
completely  local,  and  the  streaming  step  is  uniform  and  requires  little  computational  effort.  Equation  (1.7) 
is  explicit,  easy  to  implement,  and  straightforward  to  parallelize. 


1.2.  Boundary  Condition  on  a  Solid  Surface.  To  date,  most  Neumann  type  boundary  conditions 
for  a  solid  boundary  used  in  the  LBE  method  are  based  upon  the  bounce-back  boundary  condition:  A 
particle  colliding  with  a  stationary  wall  simply  reverses  its  momentum.  Much  of  the  previous  work  on  LBE 
boundary  conditions  is  devoted  to  the  analysis  and  improvement  of  the  bounce-back  boundary  condition 
[29,  10,  17,  2,  22,  7,  11,  15,  19,  30].  The  bounce-back  boundary  condition  can  attain  second-order  accuracy 
if  the  boundary  is  fictitiously  placed  half-way  between  two  nodes.  That  is,  the  second-order  accuracy  of  the 
bounce-back  boundary  condition  can  only  be  achieved  when  the  boundaries  are  located  right  in  the  middle 
of  two  neighboring  lattices  [A  =  0.5;  see  Eq.  (1.8)].  (Readers  are  referred  to  our  recent  work  [20]  for  a 
summary  of  the  previous  work.)  This  prevents  the  direct  application  of  the  bounce-back  type  boundary 
conditions  to  simulate  a  solid  body  with  smooth  curvature.  To  circumvent  this  difficulty,  Mei  &  Shyy  solved 
Eq.  (1.2)  in  curvilinear  coordinates  using  a  finite  difference  method  to  solve  for  fa  [21].  One  can  also  use 
body-fitted  curvilinear  coordinates  with  interpolation  throughout  the  entire  mesh,  except  at  the  boundaries 
where  the  bounce-back  boundary  condition  is  used  [12].  In  more  recent  works  [8,  20],  Cartesian  coordinates 
are  adopted  with  interpolation  used  only  at  the  boundaries.  These  techniques  rely  on  the  freedom  of  using 
interpolation  techniques.  We  used  the  latter  technique  in  the  present  work. 

As  shown  in  Fig.  1  for  a  2-D  projection  involving  a  3-D  body,  the  streaming  step  requires  the  knowledge 
of  fa(xb,  t),  in  which  e5  =  -ea,  at  Xb  on  the  solid  side  in  order  to  compute  fa(xf,  t)  for  the  lattice  node 
located  on  the  fluid  side  at  Xf  =  Xb  +  eaSt.  Defining 


(1.8) 


^  _  ||*/  xw 

II*/  -Xb\ 


as  the  fraction  of  an  intersected  link  in  the  fluid  region,  it  is  seen  that  0  <  A  <  1  and  the  horizontal  or 
vertical  distance  between  Xb  and  xw  is  (1  —  A )Sx  on  the  cubic  lattice. 

Based  on  the  work  of  Filippova  and  Hanel  [8],  hereinafter  referred  to  as  FH,  Mei  et  dl.  [20]  proposed 
the  following  treatment  for  on  curved  boundaries, 


(1.9) 
with 

(1.10) 
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fs(xb,  t)-( l-  x)fa(xf,  t)  +  xfa](xb,  t)  -I-  2w„p— ea  ■  uw 


faHxb,  t )  =  Wap 


,  3  9  ,  ,2  3  .  . 

If  c2  ea‘Ubf  +  2c4  (ea'Uf)  2tj2  («/■«/) 
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p - " 

Xf 
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Fig.  1.  2-D  projection  of  the  layout  of  the  regularly  spaced  lattices  and  curved  wall  boundary.  The  thick  curve  marks  the 
boundary  location.  The  solid  circles  (•)  mark  the  positions  where  particle-boundary  collision  occurs.  The  empty  (o)  and  shaded 
(•)  circles  are  fluid  sites  and  solid  sides,  respectively. 


and 

(1.11a) 

(A  -  1)  1 

Ubf  ~  ^  Uf  +  ^ Uw 

and  x  —  —(2  A  —  1) 

T 

for  A  >  1/2, 

(1.11b) 

Ubf  =  Uff  =  Uf(xf  +  eaSt,  t) 

,  (2A  -  1) 

“d  *=  (r  —  2) 

for  A  <  1/2. 

It  is  noted  that  Eq.  (1.11b)  for  Ubf  and  x  differs  from  that  originally  proposed  by  FH.  The  choice  for  uyf 
given  by  Eq.  (1.11b)  improves  the  computational  stability  for  r  <  1  and  A  <  1/2  [20].  Since  Eqs.  (1.9)  - 
(1.11)  are  in  vector  form,  they  can  be  directly  extended  to  3-D  flows  with  curved  boundaries. 

1.3.  Scope  of  the  Present  Work.  The  present  study  examines  two  issues  in  3-D  incompressible  fluid 
dynamics  simulations  with  arbitrary  boundaries  using  the  LBE  method:  i)  The  performance  of  various  3-D 
athermal  LBE  models  for  viscous  flows,  and  ii)  the  efficacy  and  reliability  of  the  extension  of  the  curved 
boundary  treatment  from  2-D  to  3-D  flows.  We  focus  on  the  stability  and  accuracy  of  the  computation  and 
the  robustness  in  handling  an  arbitrary  curved  geometry.  In  Section  2,  a  modification  of  the  choice  of  Ubf 
and  the  expression  for  x  when  A  >  1/2  is  proposed  in  order  to  further  improve  the  computational  stability 
of  the  boundary  treatment.  In  Section  3,  numerical  results  for  four  3-D  steady  flows  are  examined  and 
various  computational  issues  are  addressed.  These  four  cases  are:  (i)  pressure  driven  fully  developed  flow  in 
a  square  duct;  (ii)  3-D  lid-driven  cavity  flow;  (iii)  pressure  driven  fully  developed  flow  in  a  circular  pipe;  and 
(iv)  uniform  flow  over  a  sphere.  In  cases  (i)  and  (iii),  the  LBE-based  numerical  solutions  can  be  compared 
with  known  exact  solutions  so  that  the  accuracy  of  the  LBE  solutions  can  be  determined.  The  difference  in 
these  two  cases  is  that  A  is  a  constant  in  the  square  duct  while  A  varies  around  the  solid  boundary  in  the 
circular  pipe.  In  the  lid-driven  cavity  flow,  the  singularity  at  the  corners  between  the  moving  and  stationary 
walls  allows  for  a  performance  assessment  of  various  LBE  schemes.  The  flow  past  a  sphere  is  an  external 
flow  around  a  3-D  blunt  body.  In  all  four  cases,  detailed  assessments  are  made  in  terms  of  error  norms  and 
velocity  profiles.  It  will  be  demonstrated  that  accurate  and  robust  solutions  are  obtained  using  the  newly 
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proposed  boundary  conditions  along  with  the  selected  LBE  models. 


2.  Modification  of  the  Boundary  Condition  for  A  >  1/2.  Equations  (1.9)  -  (1.11)  are  first  applied 
to  a  fully  developed  pressure  driven  2-D  channel  flow  by  using  the  3-D  LBE  model  D3Q19.  At  the  inlet 
(i  =  1)  and  exit  (i  =  Nx,  in  which  Nx  is  the  number  of  lattices  in  the  x-direction)  the  following  zero  derivative 
condition  is  imposed  after  the  collision  step, 

(2.1a)  fa(i  =  1,  j,  k)  =  fa(i  =  2,  j,  k ), 

(2.1b)  fa(i  =  Nx,  j ,  k)  =  /„(*  =  JV*  -  1,  j,  k). 

At  k  =  1  and  k  =  Nz,  the  same  is  imposed, 

(2.2a)  fa(i,  j ,  k  =  1)  =  fa(i,  j ,  k  =  2), 

(2.2b)  fa(i ,  j,  k  =  Nz)  =  fa(i ,  j,  k  =  Nz-  1). 


The  constant  pressure  gradient  Vp  along  the  ^-direction  is  treated  as  a  body  force  and  is  included  in  the 
solution  procedure  after  the  collision  step  and  the  enforcement  of  the  above  zero-derivative  conditions  as: 


(2.3) 


t ) 


fa(Xi,  t) 


3  dp 

^  a  y  ^  a  '  X , 

c-  ax 


where  x  is  the  unit  vector  along  the  x-axis.  On  the  solid  walls  (y  =  0  and  y  =  H),  Eqs.  (1.9)  -  (1.11)  are 
used.  The  exact  solution  for  the  velocity  is  used  as  the  velocity  initial  condition.  The  equilibrium  distribution 
fiei '  function  based  on  the  exact  solution  for  the  velocity  profile  is  used  as  the  initial  condition  for  fa.  The 
pressure  gradient  is  set  to  ^  — 1.0  x  10-6.  All  computations  are  carried  out  using  double  precision. 


Fig.  2.  Stability  boundary  of  FH’s  scheme  in  a  square  duct  flow  for  A  near  1. 


Shaded  areas  under  stability  boundaries 


are  unstable  regions. 


It  was  found  that  the  computations  are  stable  for  r  close  to  0.5  (for  example,  r  =  0.505)  as  long  as  A 
is  not  too  close  to  unity  (for  example,  A  <  0.87).  When  A  is  equal  to  1,  stable  computation  can  only  be 
carried  out  for  r  no  smaller  than  0.6.  Fig.  2  shows  the  stability  bound  for  the  channel  flow  simulation  with  a 
system  size  Nx  x  Ny  x  Nz  =  5  x  35  x  5,  near  A  =  1.  Also  shown  by  the  dashed  line  is  the  stability-instability 
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boundary  for  the  channel  flow  simulation  using  D2Q9  model  and  with  a  system  size  Nx  x  Ny  =  5  x  35,  near 
A  =  1.  It  is  clear  that  similar  behavior  exists  in  both  2-D  and  3-D  channel  flow  simulations.  When  the 
computation  for  the  pressure  driven  flow  in  a  square  duct  was  carried  out  using  the  D3Q19  formulation,  a 
similar  stability  bound  was  encountered. 

Ideally,  one  would  like  to  use  a  fixed  value  of  r  for  the  entire  range  of  0  <  A  <  1  in  a  simulation. 
Computational  stability  would  then  require  the  use  of  r  around  0.6,  instead  of  a  value  that  is  close  to  0.5, 
which  makes  it  difficult  to  simulate  a  lower  viscosity,  or  higher  Reynolds  number  flow.  To  overcome  the 
restriction  imposed  by  the  numerical  stability  requirement  due  to  interpolation,  it  would  be  useful  if  one 
could  decrease  the  value  of  \  =  (2A  —  1  )/r  given  by  Eq.  (1.11a).  This  can  be  accomplished  by  using 

/  3  \  3  2A  —  1 

(2.4)  ubf  =  uf  +  —uw  and  X  =  T  +  1/2  for  A  -  V2- 

That  is,  the  velocity  ubf  is  evaluated  at  (xb  +  ea/2),  instead  of  at  xbl  using  the  information  at  Xf  and  xw 
through  linear  extrapolation. 

With  Eq.  (2.4)  replacing  Eq.  (1.11a),  the  channel  flow  simulations  using  D3Q19  lattice  model  are  carried 
out  again  for  A  from  0.85  to  1.0.  Satisfactory  results  for  the  velocity  profiles  are  obtained  for  r  =  0.505  with 
Nx  x  Ny  x  N~  =  5  x  35  x  5  in  terms  of  computational  stability.  For  A  <  0.85,  the  accuracy  of  the  solutions 
using  Eqs.  (1.11a)  and  (2.4)  is  the  same  when  the  computations  are  stable. 


3.  Results  and  Discussions. 

3.1.  Fully  Developed  Flow  in  a  Square  Duct.  For  fully  developed  flow  inside  a  square  duct  of 
height  H  defined  by  the  region  —a<y<a,  and  —a<z<a,  where  a  =  H/ 2,  the  axial  velocity  profile  can 
be  found  in  Ref.  [28,  p.  123]: 


(3.1) 


ux(y,  z) 


cosh  (««!£.) 

(2  k  +  1)« 


Figure  3  compares  the  exact  axial  velocity  profiles  at  2  =  0  and  the  LBE-based  solution  with  A  =  0.2  and 
H  =  2a  =  32.4.  A  total  of  Nx  x  Ny  x  Nz  =  13  x  35  x  35  grid  points  are  used.  The  pressure  gradient  is 
||  =  —1.0  x  10-6  and  r  =  0.52.  The  nineteen-velocity  model  is  used  in  the  simulations.  Excellent  agreement 
was  obtained. 

Figure  4a  shows  the  dependence  of  relative  Lo-norm  error, 


(3.2) 


Eo  — 


\ 


H  rH 


[wlbeO/,  z)  ~  ^exact (y,  z)]2  dydz 


0  Jo 


nH 

Kxact  (y,  z)f  dydz 


on  the  duct  height  or  the  lattice  resolution  H  =  Ny  —  3  +  2A.  The  integral  is  evaluated  by  the  trapezoidal 
rule.  As  was  demonstrated  by  Mei  et  al.  [20],  the  boundary  treatment  results  in  second  order  convergence 
for  2-D  channel  flow.  Fig.  4a  clearly  shows  that  the  total  error  (from  both  the  flow  field  and  the  boundary 
condition)  of  the  LBE  solution  in  3-D  flow  decays  quadratically. 

Figure  4b  shows  the  relative  L2-norm  error  E-2  as  a  function  of  A  in  the  duct  flow  using  13  x  35  x  35 
grid  points  and  r  =  0.52.  For  the  purpose  of  comparison,  the  relative  L2-norm  error  in  the  2-D  channel  flow 
simulation  using  the  D2Q9  model  with  Ny  =  35  and  r  =  0.52  is  also  shown.  The  relative  error  is  larger  in 
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CV2 


Fig.  3.  Comparison  of  axial  velocity  profiles  in  a  pressure- driven  square  duck  flow  at  z  =  0  between  the  exact  solution 
( solid  line)  and  the  LBE  solution  ( dashed  line)  with  A  =  0.2,  dp/dx  =  —1.0  •  10-6,  r  =  0.52,  and  H  =  32.4. 


3-D  duct  flow  than  in  the  2-D  channel  flow.  Nevertheless,  the  error  exhibits  the  same  qualitative  behavior 
in  both  2-D  and  3-D  as  a  function  of  A. 

It  should  be  noted  that  the  accuracy  of  the  D2Q9  model  and  the  D3Q19  model  is  different  in  the  sense 
that  beyond  the  conserved  moments  (density  and  momentum  in  athermal  fluids),  these  two  models  have 
different  accuracy  in  preserving  higher  order  moments  (fluxes)  [13,  14].  The  D2Q9  model  preserves  all  the 
moments  up  to  the  second  order  in  momentum  space,  which  include  momentum  fluxes,  and  maintains  the 
isotropy  of  these  moments,  whereas  the  D3Q19  model  can  preserve  density  and  momentum,  but  cannot 
maintain  the  same  accuracy  and  isotropy  of  the  fluxes  as  the  D2Q9  model.  The  only  3-D  equivalent  of  the 
D2Q9  model  in  terms  of  accuracy  of  the  moments  is  the  D3Q27  model  [13,  14], 

3.2.  Simulation  Results  for  3-D  Lid-Driven  Cavity  Flows.  Lid-driven  cavity  flow  has  been  stud¬ 
ied  extensively  in  the  CFD  community.  Most  research  has  been  focused  on  2-D  problems.  Limited  numbers 
of  reliable  numerical  results  for  steady  state  3-D  cavity  flows  have  been  obtained  in  the  past  several  years. 
In  this  study,  the  multi-block  finite  difference  solution  of  the  NS  equations  obtained  recent  by  Salom  [26]  is 
used  to  compare  with  the  present  LBE  based  results. 

The  size  of  the  cavity  is  H3 ,  the  number  of  grid  is  Nx  x  Ny  x  Nz,  and  Nx  =  Ny  =  Nz.  The  driving  lid 
is  placed  at  y  =  H,  moving  along  the  direction  of  .r-axis  with  a  speed  U  =  0.1  in  lattice  units.  Figure  5a 
compares  profiles  of  horizontal  velocity  ux{y)  obtained  using  33  x  33  x  33  lattices  with  the  solution  to  the 
NS  equations  at  x/H  =  z/H  =  0.5  for  Re  =  400.  All  three  LBE  models  (D3Q15,  D3Q19,  and  D3Q27) 
are  used.  For  the  fifteen-velocity  model,  the  computation  becomes  unstable  and  blows  up  at  this  Reynolds 
number  with  33s  lattice  resolution  and  A  =  0.5.  For  A  =  0.5,  the  nineteen- velocity  model  and  the  twenty- 
seven- velocity  model  give  very  similar  ux(y)  profiles  and  both  under-predict  slightly  the  magnitude  of  the 
minimum  in  the  profiles.  The  19-velocity  model  is  also  used  with  A  =  0.25;  there  is  a  slight  overshoot  in  the 
velocity  profiles  in  comparison  to  the  results  in  Ref.  [26].  Fig.  5b  compares  ux(y)  profiles  obtained  using  the 
fifteen-velocity  and  19-velocity  lattice  models  on  the  673  lattice  grids  and  A  =  0.5  with  the  NS  solution  [26] 
at  x/H  =  z/H  =  0.5  for  Re  =  400.  Excellent  agreement  is  observed.  Clearly,  the  nineteen- velocity  model  is 
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Fig.  4.  Error  behaviors  in  a  pres  sure- driven  square  duck  flow  with  dp/dx  =  —1.0  •  10-6,  r  =  0.52.  (a)  Dependence  of  the 
relative  L2-norm  error  on  the  lattice  resolution  H  =  Ny  —  3  +  2 A,  with  A  =  0.2.  The  straight  line  is  a  least-square  fit  of  the 
data  (symbols)  with  a  slope  of  —2.  (b)  Relative  L^-norm  error  as  a  function  of  A  for  the  3-D  duct  flow  ( solid  line)  and  the 
2-D  channel  flow  simulations. 


superior  to  the  fifteen-velocity  model.  Although  the  fifteen-velocity  model  requires  21%  less  CPU  time  and 
storage  than  the  nineteen-velocity  model  per  lattice,  it  is  not  as  robust  as  the  nineteen-velocity  model  and 
may  actually  require  more  CPU  time  and  memory  to  obtain  a  reasonable  solution  since  more  lattice  points 
are  clearly  needed. 

It  should  be  noted  that  stability  property  of  the  nineteen-velocity  model  and  the  fifteen-velocity  model 
are  significantly  different.  All  LBE  models  have  inherent  spurious  invariants  because  of  their  simple  dynamics 
[18].  However,  the  stability  of  the  LBE  models,  which  is  very  much  affected  by  these  spurious  invariants, 
differs  from  one  model  to  another,  and  also  depends  on  other  factors  such  as  boundary  conditions  and  the 
local  Reynolds  number  [18].  Among  the  three  3-D  LBE  models  (D3Q15,  D3Q19,  and  D3Q27),  the  D3Q15 
model  is  the  least  isotropic  and  therefore  is  more  prone  to  numerical  instability.  This  is  independently 
verified  in  a  recent  work  by  Kandhai  et  al.  [16].  It  was  observed  that  the  D3Q15  model  may  induce  artificial 
checkerboard  invariants  which  are  the  eigenmodes  of  the  linearized  LBGK  collision  operator  at  wave  vector 
k  =  7r;  this  can  cause  spatial  oscillations  to  develop  in  the  flow  field  at  high  Reynolds  number  [18].  Although 
it  was  pointed  out  that  the  presence  of  solid  walls  can  suppress  the  oscillation  in  certain  cases,  the  solid 
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Fig.  5.  3-D  lid-driven  cavity  flow.  Normalized  horizontal  velocity  profiles  ux(x,  y,  z)/U  at  x/H  =  z/H  =  0.5,  with 
rm  Re  =  400,  U  =  0.1.  The  Navier-Stokes  solution  (solid  line)  from  Ref.  [26]  are  compared  with  the  following  LBE  solutions. 
(a)  Nx  x  Ny  x  Nz  =  33s,  the  nineteen-velocity  LBE  solutions  with  A  =  0.5  ( dashed  line)  and  A  =  0.25  ( long  dashed  line),  and 
the  twenty-seven-velocity  LBE  solutions  with  A  =  0.5  ( dotted  line),  (b)  Nx  x  Nv  x  N~  =  673,  the  nineteen-velocity  ( dashed 
line)  and  fifteen-velocity  ( dotted  line)  LBE  solutions  with  A  =  0.5. 


walls  in  the  present  case  actually  excite  the  oscillation  by  producing  a  shear  stress  singularity  at  the  two 
corners  between  the  moving  and  stationary  walls.  Clearly,  the  D3Q19  model  is  better  suited  to  handle  flow 
singularities  than  the  D3Q15  model  in  this  case. 

Figure  6a  compares  the  profiles  of  transversal  velocity  uy(x)  obtained  from  various  3-D  LBE  models 
using  33s  lattices  (grids)  with  the  NS  solution  at  y/H  =  z/H  =  0.5  for  Re  =  400.  For  A  =  0.5  we  found 
that  the  results  from  the  twenty-seven- velocity  model  deviate  more  from  the  NS  results  of  Ref.  [26]  than 
the  results  of  the  nineteen- velocity  model  with  the  same  resolution  in  the  spatial  region  0.1  <  x/H  <  0.6. 
Both  models  under-predict  the  extrema  of  the  velocity  profile  compared  to  the  NS  solution  of  Ref.  [26].  For 
A  =  0.25,  the  results  of  the  nineteen-velocity  model  slightly  over-predict  the  extrema,  also  shown  in  Fig.  6a. 
However  the  difference  is  relatively  smaller  in  both  cases.  Figure  6b  shows  similar  results  of  velocity  profile 
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with  a  resolution  of  673  grid  points  and  the  same  Reynolds  number  Re  =  400.  With  673  lattice  resolution, 
the  result  of  the  fifteen-velocity  model  significantly  differs  from  the  result  of  the  nineteen-velocity  model  and 
that  of  the  NS  solution  in  Ref.  [26].  These  comparisons  further  suggest  that  the  nineteen- velocity  model  is 
better  than  the  fifteen-velocity  model  in  terms  of  accuracy  and  stability,  and  better  than  the  twenty-seven- 
velocity  in  terms  of  computational  efficiency.  The  nineteen-velocity  model  represents  a  good  compromise  in 
terms  of  both  computational  efficiency  and  reliability. 


x/H 


Fig.  6.  3-D  lid-driven  cavity  flow.  Normalized  transversal  velocity  profiles  uy(x,  y,  z)/U  at  y/H  =  z/H  =  0.5,  with 
rm  Re  =  400,  U  =  0.1.  The  Navier-Stokes  solution  (solid  line)  from  Ref.  [26]  are  compared  with  the  following  LBE  solutions. 
(a)  Nx  x  Ny  x  Nz  =  33s,  the  nineteen-velocity  LBE  solutions  with  A  =  0.25  ( long  dashed  line)  and  A  =  0.5  ( dotted  line),  and 
the  twenty-seven-velocity  LBE  solution  with  A  =  0.5  ( dashed  line),  (b)  Nx  x  Ny  x  Nz  =  673,  the  nineteen-velocity  ( dashed 
line)  and  fifteen-velocity  ( dotted  line)  LBE  solutions  with  A  =  0.5. 

Figures  7  and  8  show  the  effect  of  Reynolds  number  (from  100  to  2000)  on  the  profiles  of  horizontal 
velocity  ux(y)  and  transversal  velocity  uy(x)  at  x/H  =  z/H  =  0.5  based  on  the  D3Q19  model.  For  Re  =  100, 
400  and  1000,  A  =  0.5  is  used.  It  is  worth  noting  that  for  Re  =  2000,  the  system  size  of  673,  U  =  0.1,  and 
r  =  0.50325,  the  LBE  simulation  with  A  =  0.5  eventually  becomes  unstable,  although  the  steady  state  result 
of  Re  =  1000  is  used  as  the  initial  condition  for  Re  =  2000.  When  A  =  0.25  is  used  on  the  673  lattice  system, 
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no  computational  instability  occurs  and  the  steady  state  solution  is  obtained.  Weak  spatial  oscillation  in  the 
ux(y)  and  uy(x)  velocity  profiles  was  observed  for  Re  =  2000,  which  indicates  that  further  increase  in  Re 
would  require  better  spatial  resolution.  It  is  also  worth  pointing  out  that  when  FH’s  boundary  condition  [8] 
is  used  for  Re  =  2000  with  A  =  0.25,  the  solution  eventually  blows  up  even  when  converged  results  (based 
on  the  present  boundary  condition  for  A  =  0.25)  at  Re  =  2000  are  used  as  the  initial  condition. 


Fig.  7.  3-D  lid-driven  cavity  flow.  Effect  of  Reynolds  number  Re  on  the  normalized  horizontal  velocity  profiles 

ux(x,  y,  z)/U  at  x/H  =  z/H  =  0.5.  The  results  are  obtain  by  using  D3Q19  LBE  model  with  U  =  0.1,  Re  =  100, 
Nx  x  Ny  x  N~  =  33s,  A  =  0.5  ( solid  line);  U  =  0.1,  Re  =  400,  Nx  x  Ny  x  Nz  =  673,  A  =  0.5  ( dotted  line);  U  =  0.1, 
Re  =  1000,  Nx  x  Ny  x  Nz  =  673,  A  =  0.5  ( dashed  line);  and  U  =  0.1,  Re  =  2000,  Nx  x  Ny  x  Nz  =  673,  A  =  0.25  ( long  dashed 
line) . 


x/H 

Fig.  8.  3-D  lid-driven  cavity  flow.  Effect  of  Reynolds  number  Re  on  the  normalized  transversal  velocity  profiles 

uy(x,  y,  z)/U  at  y/H  =  z/H  =  0.5.  The  results  are  obtain  by  using  D3Q19  LBE  model  with  U  =  0.1,  Re  =  100, 
Nx  x  Ny  x  Nz  =  33s,  A  =  0.5  ( solid  line);  U  =  0.1,  Re  =  400,  Nx  x  Ny  x  N~  =  673,  A  =  0.5  ( dotted  line);  U  =  0.1, 
Re  =  1000,  Nx  x  Ny  x  Nz  =  673,  A  =  0.5  ( dashed  line);  and  U  =  0.1,  Re  =  2000,  Nx  x  Ny  x  Nz  =  673,  A  =  0.5  ( long  dashed 
line) . 
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3.3.  Fully  Developed  Flows  inside  a  Circular  Pipe.  Figure  9  shows  the  discretized  domain  and 
the  boundary  nodes  Xb  (denoted  by  solid  symbols)  for  flow  inside  a  circular  pipe  of  radius  R  =  9.5  lattice 
units.  Geometrically,  the  LBE  simulation  of  the  pipe  flow  differs  from  that  of  the  duct  flow  in  that  the 
fraction  of  the  intersected  link  A  is  not  constant  over  the  entire  boundary.  As  seen  in  Fig.  4b,  computational 
error  can  vary  with  A  in  the  duct  flow  and  the  difference  in  the  error  can  easily  be  as  large  as  a  factor  of  four 
for  0  <  A  <  1.  Furthermore,  the  error  is  the  smallest  when  A  is  between  0.3  to  0.6.  Hence,  it  is  reasonable 
to  expect  that  the  overall  error  in  the  solution  will  depend  on  the  distribution  of  A  in  the  entire  set  of  A. 


Fig.  9.  Schematic  for  the  boundary  nodes  x f,  ( solid  symbols)  in  yz  plane  in  the  pres  sure- driven  in  a  pipe  of  radius  R  =  9.5. 


Figure  10  shows  the  relative  Lo-norm  error  for  the  axial  velocity  profile  defined  as 


(3.3) 
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where  f l  is  the  set  of  the  discrete  lattice  grids  inside  the  pipe,  as  a  function  of  radius  R  for  R  =  3.5,  4.5, 
5.5,  9.5,  13.5,  18.5  and  23.5.  The  pressure  gradient  is  =  -1.0  x  10-6  and  r  =  0.52.  It  is  noted  that 
each  simple  summation  in  Eq.  (3.3)  is  slightly  less  than  the  exact  integration  over  the  entire  circle  due  to 
the  discretization.  To  ensure  that  such  a  treatment  does  not  affect  the  qualitative  behavior  of  the  error 
measurement,  the  centerline  axial  velocity,  uc,  is  also  compared  with  the  exact  solution  and  the  error  is 
defined  as: 


(3.4) 


Er  = 


l^c,  LBE  exact) 

exact  | 


It  is  seen  that  Ec  behaves  very  similarly  to  E2  and  both  are  non-monotonic.  This  oscillatory  behavior  could 
be  due  to  the  difference  in  the  distribution  of  A,  which  in  turn  results  in  the  difference  of  the  dissipation 
due  to  the  interpolation  around  the  boundary.  Shown  also  in  Fig.  10  is  the  error  E2  of  the  square  duct  flow 
solution  (with  A  =  0.2)  as  a  function  of  equivalent  radius  H/ir1/2,  which  exhibits  a  quadratic  convergence. 
Despite  the  non-monotonic  behavior,  it  can  still  be  seen  that  on  average,  E2  and  Ec  decay  quadratically 
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with  increasing  radius  and  the  accuracy  in  the  pipe  flow  simulation  is  comparable  to  that  in  the  square  duct 
flow  simulation. 


Fig.  10.  Variation  of  relative  errors  E2  and  Ec  in  velocity  profiles  as  a  function  of  the  pipe  radius  R  or  the  equivalent 
radius  H/ir1/2  for  a  square  duct.  Shown  are  the  Bi(R)  (n)  and  EC(R)  (O)  for  the  pres  sure- driven  pipe  flow,  and  E2(H/n 1//2) 
(-{■)  for  the  pressure- driven  square-duct  flow.  The  straight  line  is  a  least-square  fit  of  E2(H / 7T1/2)  with  a  slope  of  —2. 

Figure  11  shows  the  axial  velocity  profiles  in  the  pipe  for  R  =  3.5,  5.5,  9.5,  and  13.5  in  comparison  with 
the  exact  solution.  Even  for  a  very  small  radius  R  =  3.5,  the  LBE  solution  agrees  with  the  exact  solution 
remarkably  well.  A  noticeable  discrepancy  in  the  velocity  profile  at  R  =  9.5  is  also  observed  in  E2  and  Ec 
shown  in  Fig.  10. 


Fig.  11.  Comparison  of  the  axial  velocity  profiles  between  the  LBE  solutions  ( dashed  lines)  and  the  analytic  solutions 
(solid  lines)  for  the  pressure- driven  pipe  flow.  Shown  in  the  figure  ( from  left  to  right)  are  solutions  with  R  =  3.5,  5.5,  9.5,  and 
13.5.  dp/dx  =  -1.0  •  10“6,  r  =  0.52 


3.4.  Simulation  Results  for  a  Uniform  Flow  over  a  Sphere.  The  conventional  LBE  scheme  uses 
uniform  meshes.  Without  local  mesh  refinement,  it  is  difficult  to  compute  the  external  flow  over  a  blunt 
body  efficiently  since  a  large  number  of  grid  points  in  the  far  field  will  be  wasted.  As  a  first  attempt,  the 
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flow  over  a  sphere  is  computed  within  a  finite  region  in  the  transversal  directions. 

As  shown  in  Fig.  12,  the  outer  boundary  is  placed  at  y  =  ±H/2  and  z  =  ±H / 2.  At  y  =  —H/2,  the 
lattice  is  j  =  2.  The  boundary  conditions  at  j  =  1  for  fa’a  are  given  by  the  following  linear  extrapolation: 

(3.5)  fa(i,  1,  k)  =  2 fa(i,  2,  k)  -  fa(i,  3,  k). 

The  velocity  at  j  =  2  is  set  as 

(3.6)  u(i,  2,  k)  =  u(i ,  3,  k). 

Similar  treatment  is  applied  at  y  =  +H/2  and  z  =  ±17/2.  The  extrapolation  condition  given  by  Eqs.  (3.5) 
and  (3.6)  allow  the  flow  to  leave  the  outer  boundary.  This  helps  to  reduce  the  effect  of  the  outer  boundary 
on  the  flow  field  and  on  the  drag  force.  At  the  inlet,  a  uniform  velocity  profile  is  imposed  at  i  =  1.5  (half  way 
between  the  first  and  second  lattice  points)  and  Eq.  (1.9)  is  applied  to  obtain  the  condition  for  fa(  1,  j,  k ) 
with  x  —  0-  At  the  exit,  a  simple  extrapolation  is  used, 

(3.7)  fa(Nx,  j ,  k)  =  2 fa(Nx  -  1,  j,  k)  -  fa(Nx  -  2,  j ,  k). 

On  the  surface  of  the  sphere,  Eqs.  (1.9),  (1.10),  (1.11b),  and  (2.4)  proposed  in  this  work  are  used  to  update 
the  boundary  conditions  for  fa’s.  Only  the  LBE  nineteen- velocity  model  is  used  to  simulate  the  flow  over  a 
sphere. 


Fig.  12.  Schematic  for  computational  domain  in  the  uniform  flow  past  a  sphere. 


Figure  13  shows  the  velocity  profile  ux{y)  based  on  a  series  of  computations  carried  out  for  several  values 
of  the  radius  R  =  3.0,  3.2,  3.4,  3.6,  3.8  and  4.0  for  H/R  =  10  at  Re  =  10.  The  results  are  obtained  with 
r  =  0.7.  Fig.  14  compares  the  axial  velocity  profile  (at  y  =  z  =  0)  for  the  same  set  of  parameters.  It  is 
worth  noting  that  the  present  LBE  computation  does  not  have  sufficient  resolution  for  the  given  Reynolds 
number.  Yet  the  velocity  profiles  agree  with  each  other  accurately.  The  fact  that  we  have  obtained  a  spatially 
accurate  solution  over  a  range  of  radii  strongly  suggests  that  the  present  boundary  condition  treatment  for 
curved  geometry  in  the  LBE  method  is  capable  of  handling  more  complex  geometries  while  maintaining  good 
accuracy. 

4.  Concluding  Remarks.  Three  3-D  LBE  models,  including  the  fifteen- velocity,  the  nineteen-velocity, 
and  the  twenty-seven- velocity  model,  have  been  assessed  in  terms  of  efficiency,  accuracy,  and  robustness  in 
lid-driven  cavity  flow.  While  accurate  3-D  results  can  be  obtained  by  using  various  LBE  models,  the  nineteen- 
velocity  model  is  found  to  be  the  best  for  the  cases  investigated.  The  fifteen-velocity  model  exhibits  velocity 
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Fig.  13.  Uniform  flow  over  a  sphere  at  Re  =  10,  r  =  0.7,  and  H/R  =  10.  Comparison  of  the  normalized  velocity  profiles 
Ux(y)/U  at  x  =  z  =  0  for  various  values  for  the  sphere  radius  R  =  3.0  (x),  3.2  (A),  3.4  (+),  3.6  (O),  3.8  (dashed  line),  and 
4.0  (solid  line). 


Fig.  14.  Uniform  flow  over  a  sphere  at  Re  =  10,  r  =  0.7,  and  H/R  =  10.  Comparison  of  the  normalized  centerline  at 
velocity  profiles  ux(x)/U  (at  y  =  z  —  0)  for  various  values  for  the  sphere  radius  R  =  3.0  (dash-dot  line),  3.2  (dotted  line),  3.4 
(dashed  line),  3.6  ( dash  dot  dot  dot  line),  3.8  ( long  dashed  line),  and  4.0  (solid  line). 


oscillations  and  is  prone  to  computational  instability.  The  more  complicated  twenty-seven-velocity  model 
does  not  necessarily  give  more  accurate  results  than  the  nineteen-velocity  model  with  the  same  spatial 
resolution.  In  this  study,  we  have  also  modified  the  boundary  condition  treatment  for  LBE  method  proposed 
by  Filippova  &  Hanel  [8]  and  Mei  et  al.  [20]  when  the  fraction  of  the  intersected  link  on  the  boundary  A  is 
greater  than  one  half.  This  improves  the  computational  stability  when  A  is  close  to  1  and  r  close  to  1/2. 

The  simulations  for  flows  in  a  square  duct  and  in  a  circular  pipe  indicate  that  the  current  boundary 
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condition  treatment  for  curved  geometries  results  in  second-order  accuracy  in  3-D  flows.  The  velocity  profiles 
for  flow  over  a  sphere  show  good  self-consistency  of  the  solution  over  a  range  of  sphere  radii  used. 
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Appendix  A.  3-D  LBE  Models. 

The  D3Q15  model  has  the  following  set  of  discrete  velocities: 


'  (0,0,0), 

a 

=  0, 

rest  particle, 

(A.l)  ea 

=  < 

(±1,  0,  0)c,  (0,  ±1,  0)c, 

(0,  0,  ±l)c,  a 

=  1,  2,...,  6, 

group  I, 

l  (±1,  ±1,  ±1) 

c, 

a 

=  7,  8,...,  14, 

group  III, 

and  the  weighting  factor  wa  is  [24]: 

'  2/9, 

a  =  0, 

rest  particle, 

(A.2) 

Wa  = 

<  1/9, 

a  =  1,  2,...,  6, 

group  I, 

.  1/72, 

a  =  7,  8,...,  14, 

group  III. 

The  D3Q19  model  has  the  following  set  of  discrete  velocities: 

'  (0,0,0), 

a  =  0, 

rest  particle 

(A. 3)  ea  =  < 

(±1,  0,  0)c,  (0,  ±1,  0)c,  (0,  0,  ±l)c, 

a  =  l,  2,...,  6, 

group  I, 

.  (±1,  ±1,  0 )c,  (0 

±1,  ±l)c, 

(±1,  0,  ±l)c, 

00 

T— 1 

oo" 

I- 

II 

group  II, 

and  the  weighting  factor  wa  is  [14]: 

'  1/3, 

a  =  0, 

rest  particle, 

(A.4) 

Wa  = 

1/18, 

ol  =  1,  2,...,  6, 

group  I, 

.  1/36, 

00 
i— 1 

go" 

II 

S 

group  II. 

The  D3Q27  model  has  the  following  discrete  velocities: 

'  (0,0,0), 

a  =  0, 

rest  particle 

(A. 5)  ea  =  < 

(±1,  0,  0)c,  (0,  ±1,  0)c,  (0,  0,  ±l)c, 

a=l,  2,...,  6, 

group  I, 

(±1,  ±1,  0 )c,  (0 

±1,  ±l)c, 

(±1,  0,  ±l)c, 

00 
T— 1 

oo" 

II 

S 

group  II, 

(±1,  ±1,  ±l)c, 

a  =  19,  20,..., 

26,  group  III, 

and  the  weighting  factor  wa  is  [14]: 

8/27, 

a  =  0, 

rest  particle, 

(A. 6) 

Wa  =  < 

2/27, 

a  =  1,  2,...,  6, 

group  I, 

1/54, 

a  =  7,  8, . . . ,  18, 

group  II, 

1/216, 

a  =  19,  20, ... ,  26,  group  III  . 

In  the  above,  c  =  Sx/St ,  Sx  and  St  are  the  lattice  constant  and  the  time  step  size,  respectively.  The  lattice 
structures  for  the  D3Q15  and  D3Q19  models  are  shown  in  Fig.  15. 
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Fig.  15.  Discrete  velocity  sets  {ea}  for  the  D3Q15  (left)  and  D3Q19  (right)  models. 
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